package net.yegong.matrix;

import org.netlib.lapack.Sgesdd;
import org.netlib.util.intW;

/**
 * @author cooper
 */
public class SingularValueDecomposition {

	public SingularValueDecomposition(Matrix mat) {
		this.mat = mat;
	}

	public void calc() {
		int m = mat.nRow;
		int n = mat.nCol;
		int min = Math.min(m, n);
		float[] a = MatrixHelper.copy(mat.array);
		float[] ss = new float[min];
		float[] uu = new float[m * min];
		float[] vt = new float[min * n];
		float[] work = new float[1];
		int[] iwork = new int[8 * min];
		intW info = new intW(0);
		Sgesdd.sgesdd("S", m, n, a, 0, m, ss, 0, uu, 0, m, vt, 0, min, work, 0, -1, iwork, 0, info);

		int lwork = (int) (work[0] + 0.1);
		work = new float[lwork];

		Sgesdd.sgesdd("S", m, n, a, 0, m, ss, 0, uu, 0, m, vt, 0, min, work, 0, lwork, iwork, 0, info);
		s = new DiagonalMatrix(ss);
		u = new Matrix(uu, m, min, Matrix.FORTRAN);
		v = new Matrix(vt, min, n, Matrix.FORTRAN);
		v.transposeMe();
	}

	public void calcS() {
		int m = mat.nRow;
		int n = mat.nCol;
		int min = Math.min(m, n);
		float[] a = MatrixHelper.copy(mat.array);
		float[] ss = new float[min];
		float[] uu = null;
		float[] vt = null;
		float[] work = new float[1];
		int[] iwork = new int[8 * min];
		intW info = new intW(0);
		Sgesdd.sgesdd("N", m, n, a, 0, m, ss, 0, uu, 0, m, vt, 0, min, work, 0, -1, iwork, 0, info);

		int lwork = (int) (work[0] + 0.1);
		work = new float[lwork];

		Sgesdd.sgesdd("N", m, n, a, 0, m, ss, 0, uu, 0, m, vt, 0, min, work, 0, lwork, iwork, 0, info);
		s = new DiagonalMatrix(ss);
	}

	/**
	 * @return 2-norm
	 */
	public float norm() {
		if (s == null)
			calcS();
		return s.diag[0];
	}

	/**
	 * @return Matrix cond from 2-norm
	 */
	public float cond() {
		if (s == null)
			calcS();
		float[] diag = s.diag;
		float max = diag[0];
		float min = diag[diag.length - 1];
		return max / min;
	}

	public Matrix getU() {
		if (u == null)
			calc();
		return u;
	}

	public Matrix getV() {
		if (v == null)
			calc();
		return v;
	}

	public DiagonalMatrix getS() {
		if (s == null)
			calc();
		return s;
	}

	protected Matrix mat;
	protected Matrix u;
	protected DiagonalMatrix s;
	protected Matrix v;
}
